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^ ■ Abstract 

Q ■ The aim of this paper is to propose a 2D computational algorithm for model- 

ing of the trigger and the propagation of shallow landslides caused by rainfall. 
^ . We used a Molecular Dynamics (MD) inspired model, similar to discrete el- 

' ement method (DEM), that is suitable to model granular material and to 

observe the trajectory of single particle, so to identify its dynamical proper- 
ties. We consider that the triggering of shallow landslides is caused by the 
decrease of the static friction along the sliding surface due to water infiltra- 
tion by rainfall. Thence the triggering is caused by two following conditions: 
(a) a threshold speed of the particles and (b) a condition on the static fric- 
tion, between particles and slope surface, based on the Mohr-Coulomb failure 
criterion. The latter static condition is used in the geotechnical model to es- 
timate the possibility of landslide triggering. Finally the interaction force 
between particles is defined trough a potential that, in the absence of ex- 
perimental data, we have modeled as the Lennard- Jones 2-1 potential. In 
the model the viscosity is also introduced and for a large range of values of 
the model's parameters, we observe a characteristic velocity pattern, with 
acceleration increments, typical of real landslides. The results of simulations 
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are quite promising: the energy and the time triggering distributions of lo- 
cal avalanches shows a power law distribution, analogous to the observed 
Gutenberg-Richter and Omori power law distributions for earthquakes. Fi- 
nally it is possible to apply the method of the inverse surface displacement 
velocity [4] for predicting the failure time. 



1. Introduction 

The prediction of landslides, in particular the discovering of the triggering 
mechanism, is one of the challenging problems in earth science. The term 
landslide has been defined in the literature as a movement of a mass of rock, 
debris or earth down a slope under the force of gravity . Landslides 

occur in nature in very different ways. It is possible to classify them on the 
basis of the involved material and the type of movement [l8^ . Landslides can 
be triggered by different factors but in most cases the trigger is an intense or 
long rain. Rainfall-induced landslides involve different fields, such as engi- 
neering geology, soil mechanics, hydrology and geomorphology |]|. With the 
rapid development of computers and advanced numerical methods, detailed 
mathematical models are increasingly being applied to the investigation of 
complex process dynamics such as fiow-like landslides or debris fiows. In the 
literature, two approaches have been proposed to evaluate the dependence 
of landslides on rainfall measurements. The first approach relies on dynam- 
ical models while the second is based on the definition of empirical rainfall 
i;hresholds above which the triggering of one or more landslides is possible 



Several methods have been developed to simulate the propagation 



of a landslide; most of the numerical methods are based on a continuum 
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approach using an Eulerian point of view |2|, llO| . An alternative to these ap- 
proaches is to use Lagrangian discrete-particle methods which represent the 
material as an ensemble of interacting elements, called particles or grains. 
The commonly adopted term for the numerical methods for discrete systems 
made of non deformable elements, is the discrete element method (DEM) 
and it is particularly suitable to model granular materials, debris ows and 
ow-like landslide The DEM is very closely related to molecular dynamics 
(MD), the former method is generally distinguished by the inclusion of ro- 
tational degrees-of-freedom as well as stateful contact and often complicated 
geometries. As usual, the computational load can be very onerous with the 
increasing of the complexity or of the number of the individual element. The 
inclusion of a more detailed description of the units allows for more realistic 
simulations. However, the accuracy of the simulation has to be compared 
with the experimental data available. While for laboratory experiments it is 
possible to collect very accurate data, this is not possible for real landslides. 
These arguments motivated us in exploring the consequences of reducing the 
complexity of the model as much as possible. In this paper we present a 
toy model applied to the study of the starting and progression of particles 
down a slope, whose displacement is induced by a rainfall [9]. The inclusion 
of the effect of fluids on a granular material is a challenging problem. The 
main hypothesis of our model is that the static friction decreases as a result 
of the rain, which acts as a lubricant: this friction law is inspired by Jop 
et al. in 2006 [6]. At present we do not pretend to be able to simulate a 
real landslide or debris flow, rather we want to explore a new alternative 
approach useful for this kind of problems. The resulting numerical method. 
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similar to that of molecular dynamics (MD), is based on the use of an interac- 
tion potential, i.e. the 2-1 Lennard- Jones one. This approach is particularly 
suited for the inclusion of nonlinear terms such as those given by instanta- 
neous change of velocities, constitutive relations among different quantities, 
chemical reactions, etc. This flexibility was also exploited in the modeling 
of continuous material by means of mesoscale models. Although the model 
is still schematic, and known constitutive relations are not yet included, its 
emerging behavior is quite promising. The results are consistent with the 
behavior of real shallow landslides induced by rainfall. Emerging phenomena 
such as fractures, detachments and arching can be observed in the simula- 
tions. In particular, the model reproduces well the time distribution of local 
avalanches into the landslide, analogous to the observed Omori distributions 
for earthquakes. These power laws are in general considered the signature 
of self-organizing phenomena. As in other models, this self-organization is 
related to a large separation of time scales. The main advantage of these 
particle methods is given by the capability of following the trajectory of a 
single particle, possibly identifying its dynamical properties. 

2. The Model 

We are interested in modeling superflcial landslides, therefore we describe 
an inclined soil layer as a two-dimensional structure formed by a set of masses 
or blocks. The model is based on the interaction forces that act among blocks 
in the coordinates system along the surface. The triggering conditions are 
based on the modeling of Mohr- Coulomb law. The forces that act on the 
particles are: 
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Force of gravity 




gsm{a){mi + Wi{t),0) 



(1) 



where g is the gravity acceleration, a the angle of the slope (supposed con- 
stant), rrii is the dry mass, variable from block to block and Wi is the cumu- 
lative absorbed water in time, defined as: 



where (7^i{t) is the water absorbed during time because of the rainfall. 

Static Friction 

The static friction F^^"^ is given by: 

= iJ^sirrii + Wi{t))g cos{a){ij.sexp{-wot) + /J^siowi^ - exp{wot))). (3) 

The force in Eq. [3] depends on two friction terms, characterized by co- 
efficients fig and fisiow respectively the initial coefficient /x^ at t = and the 
final one jUgiow for t ^ oc, with fig > fisiow In synthesis, the effect of rain- 
fall is to decrease the friction on the sliding surface of the landslide during 
time (through the constant velocity of the exponential). Moreover the 
friction coefficients jig and fisiow vary randomly (in small increments) with 
the position, thus modeling the roughness of the sliding surface. 

Dynamic Friction 

When the block is moving, the applied force is: 

rf^ = /irf(mi + ^i(t))cos(a)(/Xrfexp(-^o*)+Mci^o^(l-exp(^oi)))-(-'i^). (4) 
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(2) 



Eq. m is similar to Eq. (|3j) , but the direction of the force direction is opposed 
to velocity. The friction coefficients (static and dynamic) are randomly as- 
signed to spatial zones, according to a Gaussian distribution. In this way we 
are modeling a rough sliding surface. Similarly, the friction coefficients /Hd 
and /iidiow vary randomly. 

Interaction forces among blocks 

The interaction force between two blocks or particles is defined trough a 
potential that, in the absence of experimental data, we model after a 2 — 1 
Lennard- Jones one (in Eq. ([5])). The justification of this choice is given in 
section simulation methodology. 
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Rij J \ Rij 



(5) 



where Rij = 1 is the equilibrium distance between two blocks, 6=1 and r is 
the distance between two blocks: r = y^(^7^^"^^P~+T^7^^^^- 

Force of cohesion 

At beginning, the system is prepared in equilibrium, that is, the blocks 
are disposed on a regular grid. We can assume, in agreement with the law 
of Mohr- Coulomb, modified by Terzaghi in 1943 [14] (see Eq. ([6])), to have 
a tension of cut, due to a cohesion force, also in a condition of zero normal 
tension: such principle is expressed as: 

Tf = c' + a tan((/)0. (6) 

In order to start the rupture, the tension of cut on the sliding surface equals 
an adhesive part c' plus a friction part a'n tan((/)'). Therefore, in analogy 
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with this method, the motion of the single block will not be initiate until 
the active forces (indicated in the third term of Eq. (6]) exceed the static 
friction threshold plus a cohesion term (that depends on the position in 
stochastic way). Moreover, we consider a speed threshold vd for the static- 
dynamic transition. Summing all up, the block is at rest under the following 
conditions: 



n 



(7) 



where F^^ and Vi represent the active forces and the speed. C[ is the term 
of cohesion, variable with the position, while jivi is the term of viscosity. For 
6=1 the potential energy dV = —Fdr becomes: 

V = -k [ (--\]dR = -k (inr + - - 1 ) , (8) 



1^ ij^'2i J \ 



in which we choose —1 as arbitrary constant of integration so to have zero 
potential energy at the equilibrium distance. 

3. Simulation Methodology 

In our simulations we consider an interaction among those particles at 
distance below a given threshold which in our units is 2^/^. At beginning 
the particles are arranged on a regular grid, i.e., at the instant t = each 
block is placed in the nodes of a regular rectangular grid and therefore every 
mass interacts with the eight blocks placed in the nearest and next-to-nearest 



(a) 



(b) 



Figure 1: (a) At t = we have an interaction among second neighbours, (b) Technique 
of ri-calcolous of interactions between particles: at each time step, the interactions are 
re-calculated for each mass within the interaction range. 

nodes (Figure la). For each time step, the interactions are re-calculated for 
each mass within the interaction range. This technique is used in molecular 
dynamics and congruent with principle of action and reaction (Figure lb). 

In our simulations and generally in MD the positions and velocities are 
updated using a first or second-order Verlet algorithm [3]. This algorithm 
allows a good numerical approximation and is very stable. It also does not 
require a large computational power as the forces are calculated once for each 
time step. When a block is in state of motion, the total force that acts on it 
is given by the sum of the active forces and the force of dynamic friction. 



We have to define a starting time of the landslide, for instance the time of 
the first block detachment. In case of uniform rainfall, it is simple to deduce 
theoretically this time, i.e., we can write, in the equilibrium conditions limit, 
for the single mass z. 




i i i 



(9) 



(10) 



F, = F^^HY,F^J-u•v,, (11) 

But since initially the masses are arranged in a regular grid, the interaction 
term is null as the term of viscosity that depends on the velocity, so: 
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(12) 

Vi = 0, 

thus 

\F^^\ = F^^ + Cl (13) 

that is: 

m*^sin(a) = mlgcos{a){iLis^w{-WoTp) + fisiowi^ - exi[){-WoTp))), 

m* = mi + AwTpn, (14) 

Tpn Tpl 

where Tp is the time of particle triggering, the number of temporal steps 
in simulation and At the amplitude of simulation step (At = 0.01). Solving 
the Eq. ( fT3l) we obtain: 



K + B 



exp{woTp), (15) 



where Aq (/x^ - fisiow)mi, Bq (/x^ - fisiow)^, A (tan(a) - iisiow)mi, 
B = (tan(a) - ^siow)^ and K = A - j£^. 

Eq. ( fT5l) is a transcendental equation solvable with numeric methods. 
An example of simulation is reported in the Figure [21 The triggering time of 
particles is variable from 80 to 180 temporal steps for a subset of particles 
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Figure 2: Triggering time of subset of particles depending on random variables (cohesion, 
friction, mass) in the coordinates system of the slope. 



depending on random variables (cohesion, friction, mass) in the coordinates 
system of the slope. In the Figure [3](a) the triggering time versus slope is 
shown for different values of cohesion. 

Actually, the sliding blocks could stop again after the first detachment, 
so our first definition of the starting time is not accurate. A more sensible 
definition for the starting time is based on the motion of the center of mass 
of the system. Since ours system is discrete, we get: 

X,{T*) - X,{T = 0) = > e. (16) 

In other words, we consider the starting time T* as the time for which the 
center of mass is displaced more than a distance from its starting position 
(assumed to be zero). See also the Figure |3](b) . 
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(a) 



(b) 



Figure 3: (a) Triggering time of particles versus slope for different values of cohesion, (b) 
Triggering time of landslides versus slope for increasing value of the threshold e (Eq. (p!6|) ). 

4. Simulation results 

In the simulations, an interesting behavior emerges by varying the coef- 
ficient of viscosity. For high values of viscosity, the model reproduces well 
the observed behavior of the slow shallow landslides, exhibiting a Gaussian 
distribution of mean kinetic energy increments (Figure [4](a)), 



where d = • Lowering the viscosity coefficient, the model exhibits a 

lognormal distribution (Figure [11(b), Eq. ( ITSl) considering x the logarithm 
of the data). With null viscosity, the simulation data can be fitted by an 
exponential (Figure [11(c), Eq. (ITSl)): 



f{x) ai exp(— d) 



(17) 



/(x) = aexp(6) 



(18) 



or a power law: 



f{x) = ax\ 



(19) 
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(c) (d) 

Figure 4: (a) Mean kinetic energy increment distribution of landslide, in case of /i = 0.05, 
shows a gaussian behavior, (b) Mean kinetic energy increment distribution, in case of 
/i = 0.01, shows a log- normal behaviour; the distribution of the logarithm of the same 
data is obviously gaussian. (c) Mean kinetic energy increment distribution, in case of 
/i = (exponential interpolation), (d) Mean kinetic energy increment distribution, in case 
of /i = (power law interpolation). 

which seems to better fit the data (Figure [4](d), Eq. ffTOl)). We measure 
also the intervals between the triggering time of local avalanches, i.e. we 
measures the time intervals between subsequent simulation steps + 1) 
for which the blocks start to move: in all cases a power law distribution is 
observed, but with the power coefficient decreasing with viscosity. This is 
consistent with the local triggering that is more frequent for observations in 
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(a) (b) (c) 

Figure 5: Power law distributions of difference of time triggering of the particles relative 
to the simulation with /i = 0.05 (a), /i = 0.01 (b) and /i = (c). 



"-^iilillp 




(a) (b) 

Figure 6: (a) An example of simulation in the inclined coordinate system: in this case 
arching phenomena have emerged (in red the still particles, in green the particles in mo- 
tion), (b) The displacements of the particle in the inclined coordinate, relative to the 
simulation reported in the Figure [6l^a). 

which values of fi are close to zero (Figure [5](a), Figure [5](b) and Figure [5](c)). 
Several authors j3, Q, Q have observed that some natural hazards such as 
landslides, earthquakes and forest exhibit a power law distribution. 

All results of distribution interpolation are reported in Tabled] and Table [2] 
using some estimator of the fitting accuracy: 
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Table 1: [Kinetic energy distribution (KDE) varying the coefficient of viscosity ja.] 



KDE 


^ = 0.05^ 


^ = 0.012 


fi = 0^ 


SSE 


2702 


1.862 




E. D. 


RD. 


R-square 


0.991 


1 


SSE 


50.4 


26.95 


Adj. R-square 


0.991 


1 


R-square 


0.947 


0.972 


RMSE 


7.583 


0.331 


Adj. R-square 


0.944 


0.970 


ai 


295.3 


370.4 


RMSE 


1.833 


1.34 


bi 


25.26 


18.07 


a 


15.61 


33.28 


Cl 


3.901 


0.515 


b 


-0.168 


-1.033 



[1] Gauss distribution of energy. [2] Log-Normal distribution of energy. [3] Distribution 
interpolation with exponential (E.D.) and power law (P.D.). 



SSE = Y^ivi - fix,))' 



i=l 

i?^ = l-|ff; SST = ±iy-y)' 

i=l 



(20) 



n — p — 1 



RMSE ' 



n — m 

where the first estimator is the Sum of Squared Residuals (SSE)^ the sec- 
ond is the Coefficient of Determination (i?^), the third is R Bar Squared 
[R^) and the last is the Root Mean Square Error (RMSE). In Figure [6](a) 
and Figure [6](b) the results of a simulation are shown. The behavior of the 
model is similar to that of real landslides: phenomena like fractures, arching 
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Table 2: [Landslide triggering time (T^) distribution varying the coefficient of viscosity ja.] 



Tr distribution ^ 


fi = 0.05 


fi = 0.01 


^ = 


SSE 


72.41 


26.31 


119.2 


R-Square 


0.9942 


0.9816 


0.9979 


Adj. R-square 


0.9937 


0.9803 


0.9976 


RMSE 


2.691 


1.324 


2.504 


a 


152 


37.06 


48.33 


b 


-2.508 


-1.485 


-1.183 



[1] Power law interpolation. 

and detachments are generated spontaneously during evolution of the system 
(Figure [6](a)). In the Figure [6](b) it is possible to observe the variations in the 
displacements of the particles. The higher displacements are observed at the 
base of the landslide, while smaller displacements and emerging phenomena, 
like arching, are observed in the bulk of the landslide. At this point it is 
possible to discuss the choice of the 2-1 Lennard- Jones potential: in our sim- 
ulations we tune the powers of potential, but the 2-1 Lennard- Jones allows 
to have a results similar to real landslide in term of velocity behavior, where 
is possible to assess the triggering, for example, with Fukuzono method j^]. 

In the Figure [71(a), Figure [71(b) and Figure [71(c) the trends of the modulus 
of the landslide mean velocity are reported. In all cases, by varying the 
coefficient of viscosity, we observe a transient with a rapid acceleration, in 
particular, for null viscosity (Figure [71(c) ), we observe the typical trend 
of rapid landslides In this case (rapid landslide), the failure time 
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(a) (b) (c) 

Figure 7: (a) Landslides mean velocity for simulation with ja = 0.05: the behaviour, after 
initial acceleration, is similar to a stick and slipe dynamics, (b) Landslides mean velocity 
for simulation with ja = 0.01: the behaviour is typical of some real cases with acceleration 
fases. (c) LandsHdes mean velocity for simulation with /i = 0, the behaviour is typical of 
some real cases with rapid acceleration fases (this behaviour is similar to rapid shallow 
landslides). 




(a) (b) 

Figure 8: (a) Inverse of mean velocity in the time simulatuion of landslide, in the square the 
interval for failure time assessment with Fukuzono method, (b) Square of the Figure [H^a): 
determination of failure time with Fukuzono method, the simulation data show a convex 
behavior. 

n 

is estimated by using Fukuzono method [4]; see the Figure [8](a) and the 
Figure [8](b). The time of triggering is calculated with the calibration of 
function: 

- = [/3(a-l)]^(t,-l)^, (21) 



16 



(a) 



(b) 



Figure 9: (a) Simulation in the coordinate system of the slope at t = 184 (in red the still 
particles, in green the particles in motion), (b) Relative mean displacement between near 
vertical layers of the particles along x axes of the slope at t = 184. 

where u is the mean velocity of landslide (i.e. the all particles in motion), 
the time of failure, t the time of simulation, while a and /3 are constant. With 
the calibration we obtain a = 0.8836, /3 = 2.6618 and U = 184. The behavior 
of this simulation is similar to real landslide jl3^. In the Figure [9](a) the status 
of landslide at t = 184 is shown, while in the Figure [9](b) the relative mean 
displacement between near vertical layers of the particles along X axes of the 
slope at i = 184 is reported. This distance is defined for particle positions 
Xjn as. 



where A^^ is the number of horizontal layer. 

Then in the Figure [10] the relative mean displacement between vertical 
particle layers Xi — X2 (initial layers) and X30 — x^i (central layers) in the 
time are shown. Note the time interval where the fracture is formed at initial 
acceleration phase and where the failure time is estimated with the inverse 
of velocity ([4]). Finally in Figure fTTT a) and Figure fTlT b) the progression of 




(22) 
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Figure 10: Relative mean displacement between vertical particle layers xl-x2 (initial layers) 
and x30-x31 (central layers) in the time: in evidence the time interval where the fracture 
is formed at initial acceleration phase and where the failure time is estimated with the 
inverse of velocity 




(a) 



(b) 



Figure 11: (a) Simulation in the coordinate system of the slope at t = 220 (in red the still 
particles, in green the particles in motion), (b) Simulation in the coordinate system of the 
slope at t = 260 (in red the still particles, in green the particles in motion) 
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simulated landslide is reported (time step t — 220 and t — 260). 
5. Conclusion 

A computational 2D mesoscopic models for shallow landslides, triggered 
by rainfall, is proposed. The latter is based on interacting particles to de- 
scribe the features of granular material along a slope, where a horizontal layer 
with thickness of one particle is arranged. For shallow instability movement 
we consider that the triggering is caused by the decrease of static friction 
along sliding surface. Particle triggering is caused by the trespassing of two 
conditions, e.g., a threshold speed of the particles and the static friction be- 
tween particles and slope surface, based on the modeling of the failure crite- 
rion of Mohr-Coulomb. For the prediction of the positions of these particles, 
after and during a rainfall, we use the Molecular Dynamic (MD) method, 
which is very suitable to simulate this type of systems. The results are quite 
satisfactory in order to claim that this type of modeling could represent a new 
method to simulate landslide triggered by rainfall. In our simulations emerg- 
ing phenomena such as fractures, detachments and arching can be observed. 
In particular, the model reproduces well the energy and time distribution of 
avalanches, analogous to the observed Gutenberg-Richter and Omori power 
low distributions for earthquakes. In particular the distribution of landslide 
mean kinetic energy shows a transition from Gaussian to power law, passing 
through lognormal to decrease the coefficient of viscosity up to zero. This 
behavior is compatible with slow landslides (high viscosity) and rapid land- 
shdes (low viscosity). The main advantage of these Lagrangian methods is 
given by the capability of following the trajectory of a single particle, possi- 
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bly identifying its dynamical properties. Finally, for a large range of model 
parmaters values, in our simulations we observed a velocity pattern, with 
acceleration increments, typical of real landslides [|l2|. 
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